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ABSTRACT 

The  description  and  comparison  of  several  explicit 
differencing  schemes  for  the  solution  of  the  inviscid,  non- 
thermal-conducting,  single-material,  one-dimensional  equations 
of  compressible  hydrodynamics  is  included  in  this  report.   The 
comparison  was  based  on  the  plots  obtained  from  a  series  of  one 
dimensional  test  problems  which  are  included  in  APPENDIX  A. 
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1.   INTRODUCTION 

A  shock  is  a  surface  discontinuity  in  a  fluid  which  moves 
through  the  fluid.   The  concept  of  a  discontinuous  change  in  the  fluid 
quantities  at  the  shock  position  is  convenient  although,  due  to  viscous 
effects,  the  discontinuity  is  actually  replaced  by  a  very  narrow  region 
in  which  the  fluid  parameters  change  rapidly.   The  presence  of  a  shock 
usually  makes  the  theoretical  treatment  of  the  fluid  behavior  quite 
difficult.   This  is  because  the  shock  is  a  moving  boundary,  along 
which  boundary  conditions  can  be  supplied,  but  whose  motion  is  not 
a  priori  known.   As  a  result,  analytical  solutions  in  the  presence 
of  variable  shocks  can  be  obtained  only  in  a  few  special  cases. 

The  numerical  solution  of  the  hydrodynamics  equations  is 
also  greatly  complicated  by  the  presence  of  shocks.  Very  smooth  flows 
can  be  numerically  handled  in  a  fairly  straight  forward  manner.   How- 
ever, with  most  simple  numerical  methods,  oscillations,  negative  energies 
and  other  difficulties  arise  in  the  presence  of  shocks.  Another  diffi- 
culty with  shocks  is  that  the  gradients  are  very  steep  and  any  refinement 
of  the  mesh  tends  to  make  the  gradients  even  steeper.   One  of  the  standard 
artifices  for  handling  shocks  is  artificial  viscosity.   This  method  has 
met  with  a  reasonable  degree  of  success.  However,  the  use  of  artificial 
viscosity  is  almost  an  art  rather  than  a  science  for  the  following 
reasons: 

1.   A  modification  of  the  original  equations  is  being  solved 
rather  than  the  original  equations  themselves. 


2.   There  are  many  different  types  of  artificial  viscosity 

e.g.  linear,  quadratic ,    combination  of  linear  and  quadratic. 
3-   The  parameters  in  all  artificial  viscosity  methods  are 

chosen  empirically. 
h.      Artificial  viscosity  is  applied  in  some  regions  and  sub- 
tracted in  others,  depending  on  the  particular  type  of 
artificial  viscosity  used. 
In  general,  the  above  make  it  very  hard  for  a  novice  to  use  codes 
which  include  artificial  viscosity. 

The  purpose  of  this  report  is  to  compare  several  differencing 
schemes  which  do  not  use  artificial  viscosity  for  solving  one  dimensional, 
fluid  flows  containing  shocks.   These  schemes  may  be  considered  logical 
extensions  of  the  familiar  Particle-in-Cell  method    to  the  case  of  a 
continuous  fluid.   The  differencing  schemes  are  compared  by  using  each 
of  them  to  approximate  the  fluid  behavior  in  some  special  cases  when 
analytical  solutions  are  known  and  then  comparing  the  results  with 
the  analytic  solutions. 

The  selected  test  problems  are  identical  to  those  used  by 

\2   3] 
Hicks     ,   who  compared  the  Lax-Wendroff  Lagrangian  scheme  with  PUFF 

(a  Lagrangian  scheme  which  uses  the  Von  Neumann-Richtmyer  technique  of 

artificial  viscosity).   Thus,  this  report  enables  a  simple  yet  thorough 

comparison  of  all  five  schemes.   The  test  problems  include  compression 

waves  and  rarefaction  waves  created  by  the  motion  of  a  piston,  the 

collision  of  two  shocks,  the  overtaking  of  one  shock  by  another  and 

the  shock  tube  problem. 


2.   PRELIMINARY  THEORY 

2.1  Lagrangian  and  Eulerian  Methods 

The  equations  of  hydrodynamics  are  usually  written  in  either 
the  Lagrangian  form  or  the  Eulerian  form.  These  two  forms  are  theoretically 
equivalent  although  truncated  finite  difference  approximations  to  them  may 
not  "be.  For  each  representation,  the  spatial  mesh  imposed  by  the  finite 
difference  scheme  divides  the  fluid  into  cells. 

In  the  Lagrangian  formulation,  the  co-ordinate  system  is  fixed 
in  the  fluid  and  moves  with  it.   The  cells  of  the  fluid  correspond  to 
elements  of  mass,  and  the  boundaries  of  the  cells  follow  the  paths  of 
the  fluid  elements.   The  equations  representing  the  behavior  of  an  inviscid 
fluid  in  a  one-dimensional  Lagrangian  coordinate  system  are: 


a^.-isu-o  (i) 


at   p   dx 
o 


H**i?-°  <3) 


and      P  =  f  (p,  I) 

where  V  is  specific  volume,  P  is  pressure,  I  is  specific  internal  energy, 

p   is  density,  p   is  initial  density,  x  is  the  Lagrangian  spatial  co-ordinate 


and  f  is  a  function.  The  first  three  equations  are  the  conservation  of 
mass,  momentum  and  energy,  respectively,  and  the  fourth  is  the  equation 
of  state. 

It  has  been  found  that  Lagrange  calculations  work  very  well 
for  certain  types  of  problems.   In  particular,  in  following  the  motion 
of  thin  plates  or  thin  ribbons  which  move  many  times  their  own  thickness, 
Lagrangian  codes  are  very  appropriate.   In  general,  fine  resolution  is 
found  in  shock  regions  at  the  expense  of  resolution  in  rarefaction  regions, 
Also,  the  motion  of  material  interfaces  is  closely  followed.  However,  two 
dimensional  Lagrangian  codes  behave  very  poorly  when  turbulence,  slip 
surfaces  and  other  severe  distortions  exist  in  the  fluid. 

The  Eulerian  formulation  describes  the  variation  with  time  of 
conditions  at  a  point  fixed  in  space.  The  fluid  can  be  thought  of  as 
flowing  through  a  mesh  of  cells  which  is  fixed  in  space.   The  equations 
in  Eulerian  form  are 


dp   dpU 

3t  +  3x"  =  ° 


oU  ^  UdU   1  dP 
3t    o^c   p  cbc 


dl  +  Udl  +  P  oU  =  Q 

oT   <5x  p  3x 


and  P  =  f(p,l) 


It  is  possible  to  transform  the  Eulerian  representation  to  the  equivalent 
Lagrangian  one  using  the  operator  —  =  ^r  +  -$—  where  —  represents  the 
total  time  variation  as  "seen"  by  a  particle  moving  with  the  fluid. 
The  important  property  of  an  Eulerian  calculation  is  that  it  proceeds 
without  difficulty  when  there  are  large  distortions  of  the  fluid. 
However,  material  interfaces  are  poorly  handled  and,  in  general,  it 
is  difficult  to  resolve  small  features  of  the  flow  moving  with  the 
fluid,  since  the  mesh  is  fixed  in  space.   The  use  of  massless  trace 
particles  can  partially  alleviate  this  difficulty.   Trace  particles 
play  no  part  in  the  hydrodynamics  of  the  problem.   Initially,  they 
are  imbedded  in  the  fluid  at  important  locations  e.g.  along  interfaces. 
The  motion  of  each  trace  particle  is  calculated  at  the  end  of  every 
time  step  using  a  locally  assigned  velocity  for  each  particle.   Thus, 
the  trace  particles  can  be  used  to  follow  the  motion  of  Lagrangian 
interfaces.   Another  disadvantage  with  Eulerian  methods  is  that  the 
mesh  size  is  not  automatically  reduced  in  areas  undergoing  compression 
as  it  is  with  Lagrangian  methods. 

The  methods  described  in  this  report  possess  features  of 
both  the  Lagrangian  and  Eulerian  representations. 

2.2  The  Particle-in-Cell  Method 

The  particle  in  cell  method  (often  abbreviated  to  PIC)  was 
developed  at  Los  Alamos  in  1955  by  Francis  Harlow  specifically  for 
the  solution  of  complicated  problems  in  fluid  dynamics.   It  is  based 
on  a  simple  physical  model  which  is  mathematically  rigorous.   PIC  is 
a  mixed  Lagrangian  and  Eulerian  method  retaining  some  of  the  advantages 


of  each.   It  competes  more  successfully  with  two  dimensional  methods 
mainly  because  it  handles  both  severe  distortions  (an  Eulerian  property) 
and  material  interfaces  (a  Lagrangian  property)  well.   A  brief  descrip- 
tion of  PIC  for  a  single  material  problem  follows. 

In  PIC,  the  region  through  which  the  fluid  moves  is  divided 
into  a  finite  mesh  of  Eulerian  cells  which  are  fixed  relative  to  the 
observer.   Associated  with  the  center  of  each  cell  is  a  velocity,  an 
internal  energy  and  the  total  mass  of  material  in  the  cell.   In  addition, 
as  shown  in  Figure  1,  the  fluid  itself  is  represented  by  Lagrangian  mass 
points,  or  particles,  which  can  move  through  the  Eulerian  mesh,  thus 
representing  the  motion  of  the  fluid.  The  density  of  a  cell  is  calcu- 
lated as  the  sum  of  the  masses  of  the  individual  particles  in  the  cell 
divided  by  the  volume  of  the  cell.   A  particle  is  assumed  to  have  the 
velocity  of  the  cell  in  which  it  resides.   The  mass  of  a  particle  does 
not  change,  once  a  problem  has  begun  and,  in  one  dimension,  all  particles 
usually  have  the  same  mass. 


•  ft 


-•— fr- 


Figure  1.   PIC  in  One  Dimension 


The  first  part  of  the  calculation  is  essentially  Lagrangian.   Equation 
(2)  is  differenced  to  yield 

(P  n  _  p  n  } 
~  n  =   n     j+1    j-1'  At 
3  J         2Ax        n 

A 

where  U.  is  the  velocity  of  the  fluid  at  time  (n+l)  At  and  at  the 

position  of  the  fluid  which  was  at  the  center  of  cell  j  at  time  n  At. 

Similarly,  I.  is  calculated  by  differencing  equation  (3).   Each  particle 
J 

is  then  displaced  or  transported  a  distance  U  At  where  U  is  an 
interpolated  value  of  fluid  velocity  at  the  Eulerian  position  of  the 
particle. 

During  the  transport,  each  particle  is  assumed  to  carry  with 
it  the  tilde  values  of  velocity  and  energy  of  the  cell  from  which  it 
came.   The  movement  of  some  of  the  particles  across  cell  boundaries 
causes  a  redistribution  of  mass,  momentum  and  internal  energy  to  take 
place.   This  procedure  is  very  much  like  a  Lagrangian  time  advance 
followed  immediately  by  a  conservative  rezone  to  the  original  undistorted 
Eulerian  grid. 

The  diffusion  of  particles  across  cell  boundaries  gives  rise 
to  a  dissipative  force  which  stabilizes  the  difference  equations  and 
gives  PIC  many  of  its  desirable  characteristics.   This  dissipative 
force  is  frequently  referred  to  as  effective  viscosity  since  it  is 
proportional  to  the  velocity  gradient.   Low  velocity  flow  has  caused 
some  difficulty  and  has  led  to  the  use  of  artificial  viscosity  in  addition 
to  the  effective  viscosity.   The  reader  should  see  section  2.U  for  a 
discussion  of  this  topic. 


8 


2.3  Shocks 

There  are  two  types  of  surface  discontinuity  in  a  fluid.   One 
is  called  a  contact  discontinuity  and  moves  with  the  fluid.  It  is  present 
at  the  boundary  "between  two  kinds  of  fluid,  and  it  may  arise  within  one 
fluid.  It  is  characterized  "by  the  continuity  of  pressure  and  of  the  normal 
velocity  component  across  it. 

The  other  type  of  discontinuity,  the  shock,  moves  relative 
to  the  fluid,  changing  the  state  of  each  element  as  it  sweeps  by.  For 
the  notation  used,  see  Figure  2  in  which  the  shock  is  moving  to  the  right. 


v 


El 

E2 

pl 

p2 

ul 

U2 

*1 

P2 

Figure  2.  A  Shock  Moving  to  the  Right 
In  general,  the  shock  speed,  v,  is  greater  than  the  material 
speed  on  either  side.  Also  U..  is  greater  than  U  due  to  the  shock  being 
compressive.   The  side  of  the  shock  front  through  which  the  gas  enters 
(the  right  hand  side  in  Figure  2)  will  be  called  the  front  side  of  the 
shock,  the  other,  the  back  side.  The  region  of  flow  behind  a  shock  front 
is  frequently  called  the  shock  wave.   The  shock  front  moves  with  supersonic 
speed  as  observed  from  the  front  side,  and  with  subsonic  speed  as  observed 


from  the  back  side.  There  are  three  equations  relating  variables  on  either 
side  of  the  shock  and  these  are  known  as  the  Rankine-Hugoniot  conditions. 
Rankine-Hugoniot  Conditions 

These  relations  are  derived  by  applying  the  laws  of  conservation 
of  mass,  momentum  and  energy  to  the  shock. 

1  m  =  p2  (v-U2)   -  Pl  (v-Ux) 

2  m  (U1-  U2)  =  P1  -  P2 

3  m  (Er  E2)  =  P^  -  P2U2 

2.k     Artificial  Viscosity 

As  stated  in  the  introduction,  any  finite  difference  scheme 
to  solve  the  hydrodynamics  equations  must  possess  some  dissipative  mechanism 
in  order  to  smooth  out  the  discontinuities  which  occur  when  shocks  are 
present.   Many  difficulties  in  the  calculation  can  arise  if  no  dissipative 
mechanism  is  present.   In  particular,  severe  oscillations  may  be  obtained. 
These  oscillations  remain  at  approximately  the  same  value  as  the  calculation 
continues  and  hence  are  bounded. 

The  first  dissipative  mechanism  was  introduced  by  Von  Neumann 

[1+] 

and  Richtmyer  in  1950     and  is  known  as  artificial  viscosity  because 

of  its  similarity  to  the  physical  phenomenon.  Essentially,  they  constructed 

an  artificial  medium  whose  final  behavior,  after  experiencing  the  shock 

waves,  was  sufficiently  like  that  of  the  original  medium  that  it  could 

be  used  in  the  hydrodynamics  equations  in  place  of  the  true  medium.   In 

their  method,  the  pressure  in  the  momentum  and  energy  equations  (2,3) 

is  increased  by  an  artificial  viscous  term  q.   In  other  words  p  is  replaced 

by  p  -h  q  where 
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and  b  is  a  constant.  This  viscous  dissipative  terra  causes  a  diffusion 
of  momentum  and  energy  -which  has  the  effect  of  changing  the  shock  from 
a  discontinuity  to  a  rapid  continuous  variation  of  the  fluid  variables 
over  a  narrow  region  a  few  mesh  spacings  wide.  This  smearing  of  the 
shock  front  dampens  oscillations  thus  enabling  calculations  to  proceed 
as  if  there  were  a  smooth  flow.   Across  the  narrow  region  occupied  by 
the  shock,  the  Rankine-Hugoniot  equations  hold.   Artificial  viscosity 
converts  kinetic  energy  into  internal  energy  across  shocks  thus  increasing 
the  entropy. 

Artificial  viscosity  works  well  in  one  dimensional  codes  but 

f5] 
has  not  worked  so  well  in  higher  dimensions.     There  are  several 

other  types  of  artificial  viscosity  in  fairly  common  use  e.g.  the 

[6] 
Landshoff  type    in  which 


q  =  q^  +  Wqi  +  (1-W)  ^ 


where  q_  is  the  Von  Neumann-Richtmyer  term 

NR 


a 
and  q  =  -  -^  pAx  U 


ll  2  ^   x 


/   %2  At  U 
^      2    p 


with  0  <  W  <  1. 
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2.5  Effective  Viscosity 

It  has  been  shewn  for  PIC  by  Taylor  expansions  of  the 
difference  equations  and  neglecting  higher  order  terms  ,  that 
solutions  to  the  PIC  finite  difference  equations  approximate 


o-U   pudu   oT   d   ,.  ,  dtK  n   N 

p  St  +   5E  +  3x"  =  35E  ix    35E>  (1+) 


and 


Si  puSi  pSu  d  ,,^l   ^,/^U\2  ,cx 


cH     oSc    3x  "  3x  ^  5x'     xix 


where     V  =  ^-Ax  PU 


more  closely  than 


n   Bu   pudu   dp 
at     ox   ox 


„  Si  pudi  pdu  _ 

md      PSt+   olT  3i=° 


Some  of  the  terms  on  the  right  hand  side  of  equations  (h)   and  (5)  are 
proportional  to  the  velocity  gradient  and  can  be  related  to  the  physical 
concept  of  viscosity.  These  additional  terms  arise  as  a  result  of  the 
mass  transport  mechanism  and  are  responsible  for  producing  the  dissipa- 
tive  effects  that  give  stability  to  the  method.  Thus,  PIC  and  its 
related  methods  are  a  class  of  codes  whose  dissipative  force  is  referred 
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to  as  effective  viscosity  or  implied  viscosity  or  smoothing.  These 
terms  allow  PIC  to  give  accurate  approximations  to  high  velocity  flow. 
The  V  terms  have  their  principal  effectiveness  in  regions  where  the 
velocity  gradient  is  large  or  where  fluctuations  occur  in  high  speed 
flow.   At  low  fluid  speeds  or  in  regions  of  near  stagnation,  the  veloc- 
ity gradients  are  small  and  hence  the  effective  viscosity  terms  are  not 
large  enough  to  make  the  dissipative  force  effective  so  that  oscillations 
develop.  These  oscillations  will  grow  until  they  become  comparable  to 
the  local  sound  speed  when  the  effective  viscosity  terms  become  large 
enough  to  prevent  any  further  growth  of  the  instabilities  (as  in  test 
problem  6b).  Thus,  the  instabilities  are  bounded.  These  difficulties 
can  be  alleviated  by  the  use  of  artificial  viscous  pressure  terms  in 
the  equations  for  momentum  and  energy  (see  section  2.k).     However, 
these  terms  should  be  added  only  where  the  fluid  speed  is  small  com- 
pared to  the  local  sound  speed.   If  applied  in  high  velocity  regions, 
the  artificial  viscosity  will  smear  out  many  features  of  interest. 

Some  authors  question  whether  the  smoothing  scheme  in  PIC 

f7  81 
can  be  directly  related  to  physical  concepts.  '  They  have  tried 

to  show  that  doing  so  results  in  many  contradictions.   For  example, 
the  values  of  the  additional  terms  far  exceed  the  values  of  the  anal- 
ogous physical  terms. 
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3.   DERIVATION  OF  THE  DIFFERENCE  EQUATIONS 

The  region  occupied  by  the  fluid  is  considered  to  be  discretized 
into  a  mesh  of  equal  sized  cells  fixed  in  space.   Each  cell  is  of  width  Ax 
and  is  assigned  some  index  j.  Half  integer  indices  are  used  to  denote 
cell  boundaries. 


Ax 


Ax 


Ax 


Figure  3.  The  Mesh  of  Cells 
The  fluid  is  then  described  at  any  instant  of  time  by  specifying  the  velocity, 
density  and  pressure  for  each  cell.   These  values  are  considered  to  be  known 
at  the  center  of  each  cell.  The  calculation  advances  explicitly  in  time 
steps  i.e.  cell  quantities  at  time  (t  +  At)  are  calculated  in  terms  of 
values  of  those  quantities  at  t. 

The  system  of  equations  to  be  solved  is 


do    5  (pU) 


£ 


5t 


=  o 


dU       UcVQ  +    idP  _ 
cH  dx       Pdx 


0 


(6) 


(7) 


lU 


||     o  *  +  i   a  (ro)    .  (8) 

at     ox   p  ax  v  y 


and  P  =  f  (P,I) 
which  represent  the  conservation  of  mass,  momentum  and  energy  and  the 
equation  of  state,  respectively,  where  the  symbols  have  their  usual  meanings, 
As  in  PIC,  the  solution  is  accomplished  by  performing  first  a  Lagrangian 
update,  then  rezoning  to  the  original  Eulerian  mesh.  During  the  Lagrangian 
phase,  the  operator  ;—  +  °r  is  replaced  by  -rr,  the  total  derivative. 
Equations  (7)  and  (8)  then  become 

dt   p  ax 


2|  +  I  |Mi  .  0  (10) 

dt   p  ox 


Equation  (10)  can  be  rewritten  as 


at  N*^  %  i  <™>   ■  ° 


ai      p  au  J  TTfaU      1  aP-,         n 

or  =7-  +  —  r—  +  UL^r  +  —  c— J      =     0 

at      p  ax         at      p  ax 


Now,    equation   (9)   implies   that 

j£+£jffl  .   0  (11) 

at      p  ax 


The  intermediate  or  tilde  values  of  velocity  and  specific  internal  energy 
at  the  end  of  the  Lagrangian  advance  are  given  by  the  following  equations. 
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3-   DERIVATION  OF  THE  DIFFERENCE  EQUATIONS 

The  region  occupied  by  the  fluid  is  considered  to  be  discretized 
into  a  mesh  of  equal  sized  cells  fixed  in  space.   Each  cell  is  of  width  Ax 
and  is  assigned  some  index  j.  Half  integer  indices  are  used  to  denote 
cell  boundaries. 


Ax 


Ax 


Ax 


-it- 
Figure  3.  The  Mesh  of  Cells 
The  fluid  is  then  described  at  any  instant  of  time  by  specifying  the  velocity, 
density  and  pressure  for  each  cell.  These  values  are  considered  to  be  known 
at  the  center  of  each  cell.  The  calculation  advances  explicitly  in  time 
steps  i.e.  cell  quantities  at  time  (t  +  At)  are  calculated  in  terms  of 
values  of  those  quantities  at  t. 

The  system  of  equations  to  be  solved  is 
So    o-  (pU) 


5t 


=  0 


(6) 


I  +  UcUJ  +  idp  = 
dx   pdx 


5t 


(7) 


lU 


||+U^+i-i(HT)   .  0  (8) 

dt     dx   p  dx  v  ' 


and  P  =  f  (P,I) 
which  represent  the  conservation  of  mass,  momentum  and  energy  and  the 
equation  of  state,  respectively,  where  the  symbols  have  their  usual  meanings, 
As  in  PIC,  the  solution  is  accomplished  by  performing  first  a  Lagrangian 
update,  then  rezoning  to  the  original  Eulerian  mesh.   During  the  Lagrangian 
phase,  the  operator  :-— ■  +  4-  is  replaced  by  -rr-,  the  total  derivative. 
Equations  (7)  and  (8)  then  become 

dt   p  dx 


M  t  1  |IM  _  (10) 

dt   p  dx 


Equation  (10)  can  be  rewritten  as 


§1 N^  ♦  ?  §1  <ro>  ■  ° 


or  |i  t  £  |H  +  VM  ,  i  |E]    .    0 

dt       p  dx  dt       p  dx 


Now,  equation  (9)  implies  that 
dt   p  dx 


The  intermediate  or  tilde  values  of  velocity  and  specific  internal  energy 
at  the  end  of  the  Lagrangian  advance  are  given  by  the  following  equations. 
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u  n  =  u.n  -  (p.nn  -  p.nJ  -^ <12> 

V  =  1°  -  Cf",-  p",)  P.i"  At           (13) 

°       '        °  3    2p.n* 

where  U.   =  —  [U.   +  U.  ] 

J      2   j  j 


These  intermediate  values  correspond  to  values  at   time  (n+l)  At  at  the 

point  in  the  fluid  that  was  at  the  center  of  cell  j  at  time  n  At. 

To  first  order,  the  tilde  quantities  are  values  at  time  (n+l)  At  at  the 

Eulerian  point  x.  +  U.  At.  The  use  of  time-centered  values  of  velocity 

—  n  ~  n 

U.   in  the  equation  for  I.  results  m  greater  stability  and  better  behavior 

of  fluid  entropy  than  if  the  simple  value  U.  were  used.     Equation  (l^) 

J 

is  used  to  calculate  the  intermediate  value  of  specific  energy. 

E.   =  I.  +  |  (U.)2  (11+) 

3  3        2   j' 


The  distorted  Lagrangian  mesh  is  now  rezoned  back  to  the  original  Eulerian 
mesh.   In  effect,  this  is  a  mass  transport  across  cell  boundaries.  The 
mass  flow  across  the  interface  between  cell  j  and  cell  (j  f  l)  is  defined 

by 

n       n    "*n 
AM.  1  =  p.  1  U.  1  At 
J+2      J+2   0+2 

*n 
The  quantity  U.  1  is  known  as  the  weighted  velocity. 

J+2 
Three  techniques  for  calculating  U.  1  define  the  three  differencing 

schemes  used  in  this  report  which  are  called  continuous  rezone,  OIL  and  donor. 

The  interface  values  of  density  and  velocity  cannot  be  calculated  by  simply 
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averaging  adjacent  cell  values,  because  the  method  obtained  is  highly- 
unstable  and  virtually  useless.   The  density  weighting  and  velocity  weighting, 
as  they  are  called,  must  be  done  with  greater  care.   The  equations  for  the 
density  weighting  and  velocity  weighting  are  simply  stated  in  this  section 
but  are  derived  at  length  in  section  3.1.  The  direction  of  the  flow  at 
the  interface  between  cell  j  and  cell  (j+l)  is  defined  by  the  sign  of 

U  1  where  U  1  =  5-  [U   +  U   ]. 
«J  p        <J    p        J      <J 

Density  Weighting  for  all  schemes 


If  U.nl  >  0      then 
J+2 

n-> 
P.  1  = 

J+2 

n 
P. 

If  U.nl  =  0    then 
J+2 

J+2 

|[p."  + 

If  U.ni  <  0      then 
^2 

PA 

Velocity  Weighting 

(a) 

continuous  rezone 

If  U, 

ni  >0 
i+5  - 

*n 
then  U.  1 

J+2 

n+1-, 

pj     ] 


u.ni 

J+2 


i  *  C^i  -  Vi>-4i 


2A  x 


^    PI 

U     1 
kL   .„         ,,  *n,  J+2 


If  U.   1  <  0         then  U.   1 


Jf2  J+2         ,        ,~  n        «  ik  At 


1  +■    (U  "     -  U/  ) 


j+2  j    '2A  x 


(b)      OIL  ^_n1 

*n.  ^+2 


U.    1 


Jfo  i   j.  fu  n        ~  nN     At 


2 


1  ♦  (u »    -  U  -) 


J+l  j    '   ^2Sc 
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density  for  the  other  two  schemes  gives  fairly  good  results. 

The  three  velocity  weighting  schemes  used  in  this  report 
are  described  below. 
3.1.1  Continuous  Rezone 
Assume  that  the  cell  in  the  diagram  is  moving  from  left  to  right. 


i.ni 


j 


Figure  k.      Continuous  Rezone 
The  undashed  rectangle  in  Figure  h     represents  the  position 
of  cell  j  at  time  n.  The  dashed  rectangle  represents  the  position  of 
the  same  rectangle  after  the  initial  Lagrangian  flow.   The  density,  specific 
internal  energy  and  momentum  are  assumed  constant  throughout  the  dashed 
cell.  The  velocity  and  specific  internal  energy  at  the  center  of  the 
dashed  cell  are  given  by  equations  (12)  and  (13)  respectively.   The  quantity 
of  mass  which  is  in  the  dashed  cell  and  has  crossed  over  into  Eulerian  cell 


(j-t-l)  during  the  rezone  is  calculated  from  the  fraction 


c '      O  _    <J 


(16) 


n 


since  the  density  over  the  dashed  cell  is  assumed  constant.  The  mass  of 
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cell  j  before  the  Lagrangian  step  equalled  M.n  and  this  obviously  equals 

the  mass  in  the  dashed  cell  after  the  Lagrangian  step.  It  can  be  seen 

that  to  first  order 

Ln     =     Ax  +  (U.nl  -  U.ni)  At 
J  J+2     D   2 

and      £nl     =     U.nl  At 
J+2      J+2 


~  n      n 
A,,  n,    U .  1  At  P  . 
Hence  AM  1  =   j+p     J 

op 

1  -H  (U*   -  U.nJ  #- 
0+1    J-l  2  Ax 

Hence  from  equation  (l6) 


2        2 


/~  n    k  n  i  At 
1  +  (U.  _  -  U.  .  )■ 
v  .1+1    .-]-!': 


j-1'2  Ax 


If  the  flow  is  from  right  to  left,  then  a  similar  analysis  yields 
*nn    U.Vu* 
J  2        2 


ns  At 


1  +  (u.nc  -  u.n)  4 

K    i+2  1  '  2 


j+2    j  '  2  Ax 


3.1.2  OIL 


The  OIL  velocity  weighting  scheme  has  been  in  use  for  some  time 
in  two  dimensional  calculations    .   The  derivation  of  OIL  is  similar 
to  that  of  continuous  rezone  except  that  the  cell  to  be  followed  is  centered 
at  the  interface  between  cells  j  and  (j+l),  and  is  Ax  in  length. 
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L.ni 
J+2 


1                                                                                   1 
1                                                                                   1 

1                1 

1 

1   \          ': 

1                        1                        1 
1                                                1 

;.  d+i    !        ;. 

1 
1 
1 

Figure  5.   OIL 
The  mass  transported  into  cell  (j+l)  by  the  Lagrangian  transport  step  is 

A..  n,   L .  1  M  . 

J+2 

L.n 


where  L.n  =  Ax  +  (U."  -  U.n)   At 


and  L.nl  =  U.nl  At 


J+- 


J+2 


*n 
Hence  U.  1 

J+2 


u.ni 

J+2 


/~  n 

1  +  (V 


u.n) 


At 
Ax 


The  same  analysis  holds  if  the  flow  is  in  either  direction.  Another  derivation 
using  Taylor  Series  expansions  can  be  found  in  reference  [10]. 


3.1.3  Donor 

The  value  of  velocity  chosen  is  simply  the  tilde  velocity  of 
the  donor  cell.   Hence 
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*n      ~  n       ~  n 
U.  1  =   U.     if  U.  1  >  0 

J+2     J         J+2 


■X-n  ~  n 

U.  1  =  0  if  U.  1  =  0 
J+2         J+2 


U.  1  =  U."  if  U.  1  <  n 
0+2     J+1     J+2 


3.1.U  Rich 

This  method      is  included  for  the  sake  of  completeness  although 
it  was  not  applied  to  the  test  problems  due  to  its  poor  stability  properties. 

The  weighted  values  of  density  and  velocity  are  obtained  from  the 
first  two  terms  of  a  Taylor  series  expansion  about  the  donor  cell.  Assuming 


that  the  flow  is  from  left  to  right  (U.nl  X)),  and  that 

•v~  ^   2     ~  n         ~  ~ 

~  n         ~  n   .    Ax     oU    l  n         ~  n      ,    U . , .,    -  U 

"2 


U.    1  =  U.      +  =~     s-      .    =     U.         +      .1  +  1  J-1  >  0 

J+O  J  2       dx     [J  J  -d r « 


then       AM.nl  =      [U.n  +  Uj  +  l"     ,j-l]      [p  .n     +  Ptj+1   "  pj-l]  At 
J+2  ^  1+  ^  1+ 


Similarly,    if  U.nl     <  0  and  if 
3  2 


-\~  ~  n         ~  n 

U^l  =  U.nn    -     %^-   |n   ,          =     U*     -  Uj+2   -  Uj        <0 
j+5-         j+1  2dx'j+l  j+1       -" j-— •L- 


~  n         ~  n  n  n 

then         AM.   1  =    [U.    .    -     j+2         j    J    [p.      -     ,]+2       Kj    J     At 

j+2         J+1         — JT~^         J  — T~^ 


Otherwise  AM.   1  =     0 
J+2 
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3.2  Conservation  of  Mass,  Momentum  and  Energy 

The  change  in  the  mass  of  cell  j  between  the  start  and  the  end 

of  each  iteration  is 

...  n   r  n+-l    n-,    . 
£M.   =  Lp  .    -  p.  Ax 
3  3  3 

=  m.ni  -    mni 

3-2  J+2 


where  the  half  integer  subscripts  indicate  quantities  moving  across  cell 

boundaries.   Therefore,  the  change  in  mass  of  the  entire  system  is 
J 

mn    =  Z  ,  [  mni  -    £M.ni] 

3=1  3 "2  J+2 

n        n 
Note  that  £M.   1  =AM//.1v   1N  because  the  mass  transported  across  the 

right  hand  side  of  cell  j  equals  the  mass  tranported  across  the  left  hand 

side  of  cell  (j+l).  Hence 

/Mn  =  ®£     -     ZM_nl 
J.       J+„ 

2         * 

Therefore,  the  only  change  in  mass  of  the  total  system  occurs  by  mass 
crossing  the  boundaries  of  the  system  and  hence,  mass  is  conserved. 

Redistribution  of  energy  and  momentum  occurs  during  the  calculation 
of  tilde  values  and  also  during  the  rezoning  from  the  Lagrangian  mesh  to  the 
Eulerian  mesh.  During  the  calculation  of  tilde  values,  the  changes  in 

energy  and  momentum  are  the  following 

J    n  r~n   _nn  A 

^  =  jLi  PJ  [  j  "  V  * 

J 

/Mom  =  Z  n      pn   [U1?     -  U1?]    Ax 
J=l        3        3  3 


2*+ 


Therefore 


~n 

AM 


=  Z   -  [P.nl  -  P*l]  At 


J-o    J '2 


As  before,  all  terms  cancel  the  first  and  the  last,  leaving 

AMn   =   (*rf  -  PTni)  At 
om       1    J+— ' 

2 


Also 


^     =     Z,      k-?*   At  (U*  -  U.nJ  +  Pn  Ax  (U*  +  Un)  (U*  -   U*)] 

t    ,i=i  2   ,i    .i+i   .i-iy   .i     .i    y   .1    .r 


J' 


=  Z,  %  [-p?  (u*  -u")  -un  (p.*  -p.nJ] 

.1=1  2    ,i   ,1+1    ,i-ly     ,i   ,i+l    .i-l 


J   J 


J   J 


At  r  n 


1   O     1   O     J   J+l     ,1   J+l 


=  At 


p?  if1  +  5?  Pn    -  p  n  if1  -  PA  u" 

1  o    1  o      o  o    J+^  J+l 

2  2 


-  uTni  p*  +  p"  u* 

J+P  J+l    J+l  J+l 


Therefore,  the  only  changes  in  momentum  and  energy  of  the  system  during 
the  calculation  of  tilde  values  occurs  by  momentum  and  energy  crossing 
the  boundaries  of  the  system.  Hence,  momentum  and  energy  are  conserved 
during  this  part  of  the  calculation. 

The  equations  for  the  rezoning  phase  are  directly  derived  from 
the  equations  for  conservation  of  momentum  and  energy  and  hence  momentum 
and  energy  must  be  conserved  during  the  rezone.  The  energy  and  momentum 
changes  of  the  system  occurring  during  the  rezone  are 
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n     ...n  ~n     A.„n  n  ~n 


omR       1  L      J+2  R 
2  d 


Thus,  the  mass,  momentum  and  energy  of  the  system  are  conserved 
identically  despite  finite  difference  approximations.  However,  corrections 
must  be  made  at  transmittive  boundaries  and  pistons. 

It  is  interesting  that  the  derivations  of  the  conservation  of 
mass,  momentum  and  energy  are  independent  of  the  velocity  weighting  scheme 
used.   This  is  because  the  conservation  laws  are  concerned  with  cell  centered 
quantities  but  the  velocity  weighting  schemes  are  defined  by  interface 
quantities. 

The  only  discrepancies  arising  between  the  theoretical  energy 
of  the  system  and  the  numerically  approximated  values  are  due  to  round- 
off errors.  The  values  of  system  energy  obtained  for  test  problem  1A 
(see  section  5.2)  which  are  typical  of  the  values  obtained  for  the  rest 
of  the  problems  are: 

Initial  system  energy  I.1363XIO  gms. 

System  energy  at  .1  seconds  minus  work  done  by  piston 

(a)  continuous  rezone       I.I36OXIO  gms. 

(b)  OIL  1.1359X109  gms. 

(c)  donor  1.1353X109  gms. 


26 


3.3  Order  of  the  Truncation  Errors 

By  applying  Taylor's  series  expansions  to  the  finite  difference 
equations  for  continuous  rezone,  OIL  and  donor,  and  comparison  with  the 
equations  of  hydrodynamics,  the  order  of  the  errors  involved  in  the  approxi- 
mations can  be  found.   The  methods  are  of  at  least  the  order  calculated  and 
might  be  higher.   The  derivations  are  carried  out  assuming  that  the  flow  is 
from  left  to  right.   Combinations  of  left  to  right  and  right  to  left  flow 
have  not  been  examined.   The  order  of  the  truncation  errors  for  all  three 
velocity  weighting  schemes  for  each  of  the  three  equations  of  hydrodynamics 

are  simply  stated  in  this  section  and  are  derived  in  detail  in  APPENDIX  B. 

2 

(a)  continuous  rezone  is  0  (a  x)  +  0   (At)  +■  0  (■—  ) 

for  each  of  the  three  conservation  equations 

2 

(b)  OIL  is  0   (Ax)  +  0   (At)  +  0  (j-   )  for  each  of  the  three 

conservation  equations 

(c)  donor  is  0   (Ax)  +0   (At)  for  the  conservation  of  mass 

At2 

and  momentum,  and  is  0   (At)  +  0   (Ax)  •*-  0  (—  )  for  the 

ZjX 

conservation  of  energy. 
Thus,  the  set  of  finite  difference  equations  representing  each 
of  the  three  schemes  continuous  rezone,  OIL  and  donor  is  of  order  0  (A  x) 
+  o   (At)  +  0  (g  ) 
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k.      TIME  STEP  CONTROL  AND  TREATMENT  OF  BOUNDARIES 

k.l     Time  Step  Control 

Three  criteria  were  used  for  calculating  the  time  step.   The 
first  was  the  familiar  Courant-Friedrichs-Lewy  stability  criterion  which 
states  that  at  any  time  step  n,  no  signal  travelling  at  the  local  sound 
speed  C  can  propagate  through  more  than  one  cell  in  one  time  step.  Hence, 
for  each  n 

Atn  1  -^ 

cn 

3 

for  all  j  in  the  region  of  interest 

.,n    „   .   rA  x-, 
or  At  =  F  min  L J 

i     c." 

< 

where  0  <  F  -  1. 

F  is  defined  as  the  Courant-Friedrichs-Lewy  number. 

The  second  limitation  which  applies  to  supersonic  flow  problems 
only,  prohibits  the  transport  of  mass  across  more  than  one  cell  in  one 
time  interval.   The  fluid  model  provides  no  mechanism  for  transfer  of 
mass  between  non-adjacent  cells.   Hence 

Atn  1  -& 

u.n 

for  all  j  in  the  region  of  interest 

The  third  limitation  is  strictly  applicable  to  the  OIL  method 
but  was  applied  to  the  other  schemes  as  well,  enabling  all  three  techniques 
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to  use  the  same  Courant-Friedrichs-Lewy  number.   Thus,  continuous  rezone 
and  donor  should  he  able  to  run  faster  than  OIL. 

Suppose  At  is  calculated  from  the  Courant-Friedrichs-Lewy  criterion. 

Then 

.,  n     <     Ax       ,  _7  /max    .    TT  n   .  _  nN 

At       —    — -    where  V  =  max   (    .  U  ,   C.    ) 

V  3      '      3  3 

FA  x  < 

or     At  =     — - —    where  F  —    1 

Assuming  that  the  flow  is  from  left  to  right,  suppose  that  some 

U.n  =  V.   Then 
3 

U.n  At  =  FAx  (28) 

3 


However,  the  fluid  model  for  OIL  prevents  the  left  hand  boundary  of  the 
cell  whose  center  is  at  (j  +  p)  from  crossing  into  cell  (j  +  l). 


Hence 


tt  n   a*.  <        A  x 

U5   At  -  —  (29) 


By  comparing  equations  (28)  and  (29)  it  can  be  seen  that  F  —  _-  is  a  necessary 

condition  in  order  to  remain  consistent  with  the  fluid  model.   It  may  also 

~  n  v.  1 

happen  that  some  U.  >  V,  in  which  case  F  <  _-  for  stability.   It  is  necessary 

to  show  that  the  condition  U.  >  V  can  exist  to  support  the  hypothesis  that 

J 

F<i. 

Suppose  that  u,    =   V  for  some  j  =  k.   It  is  known  that 


Uv  =  u?  - 
k      k 


* 


J  K 


Hence,    if  ^   |        <  0,    then  UjJ  >  V. 
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In  the  test  problems,  F  was  chosen  to  be  equal  to  .h.      This 
is  equivalent  to  assuming  that  no  U\  will  increase  by  more  than  25% 
over  V  in  one  time  step. 

Some  authors     have  experienced  difficulties  with  negative 
energies  and  negative  volumes  and  have  incorporated  additional  time 
step  controls  to  handle  them.   However,  this  was  not  found  necessary 
in  this  report. 


k.2     Treatment  of  Boundaries 

A  transmittive  boundary  at  the  interface  (j  +  -|)  is  defined 
by  the  following  equations : 


OX  '  0+2 


t —    •  1  =  u 
ox   j+i 


£  I    .1  =  0 

ox   j++ 


s'jiK0 


For  a  left  hand  boundary  at  J  =  ■§-,  the  above  equations  are  satisfied 

by  using  a  dummy  cell  with  index  0  and  letting  U  =  U  ,  P   =  p  ,  P 

P  and  I  =  I  .   A  right  hand  boundary  is  handled  in  a  similar  fashion. 
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It  is  necessary  to  calculate  the  mass  and  energy  transported  across 
a  transmittive  boundary  in  order  to  check  the  conservation  laws.   The 
conditions  at  a  reflective  boundary  are  the  same  as  above  except  U  = 
-U  .   Moving  piston  boundaries  are  complicated  by  the  work  terms  done 
on  the  gas  by  the  piston. 
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5.   DESCRIPTIONS,  NUMERICAL  VALUES  AND  DISCUSSIONS  OF  TEST  PROBLEMS 

5.1  Introduction 

The  test  problems  are  intended  to  exhibit  typical  behavior  of 

the  hydrocodes  in  a  wide  variety  of  situations  and  are  not  difficult  in 

r-5] 
themselves.   In  the  report  from  which  the  test  problems  are  taken    . 

there  are  seven  test  problems  with  a  few  variations  of  each,  which  are 
designated  A,  B,  C  etc.   The  difference  between  these  variations  is  in 
the  initial  numerical  values  of  one  or  more  of  the  fluid  parameters. 
The  test  problems  in  [3]  which  were  not  run  with  continuous  rezone,  OIL 
and  donor  are  described  in  the  next  paragraph.   The  format  for  the  descrip- 
tion of  each  test  problem  will  consist  of  a  brief  description,  initial 
numerical  values,  boundary  conditions  and  plot  times  together  with  initial 
and  final  plots  of  pressure,  density,  specific  internal  energy  and  velocity. 
Sometimes,  initial  plots  are  not  given  if  one  or  more  parameters  is  constant 
throughout  the  region  of  interest.  For  further  details  about  the  test 
problems,  the  reader  should  consult  reference  [2],  Also,  the  PUFF,  Lax 
Wendrof f  and  analytic  solutions  may  be  found  in  reference  [3  ]. 

None  of  the  last  three  parts  of  problem  two  and  none  of  problem 
four  was  run  due  to  the  presence  of  a  vacuum.   When  a  vacuum  arises,  the 
density  equals  zero  and  equations  using  density  as  a  divisor  cannot  be 
used.  The  last  part  of  problem  seven  was  not  run,  since  it  is  just  a 
scaling  down  of  the  first  part  by  a  factor  of  100  and  no  new  features 
are  introduced.  Also,  problem  three  was  not  run  mainly  because  it  is 
a  pure  compression  problem  and  in  this  respect  would  yield  little  informa- 
tion not  already  given  by  problem  one. 
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The  chief  difference  between  the  test  problems  in  Hick's  reports 
and  the  ones  in  this  report  are  that  the  latter  may  be  translated  in  the 
positive  x  direction  to  avoid  negative  or  zero  subscripts.  This  will  be 
clearly  stated  where  appropriate.  Also,  there  might  be  a  slight  difference 
between  the  total  number  of  cells  in  the  mesh  used  in  this  report  and  in 
reference  [3].   Neither  of  these  differences  will  change  the  calculated 
error  norms  because  the  cell  sizes  in  each  report  are  the  same  and  the 
discrepancy  in  the  total  number  of  cells  occurs  at  the  boundaries  of  the 
problem  where  the  exact  and  approximate  solutions  are  equal. 
Errors 

The  analytic  solution  to  the  mathematical  model  of  the  physical 
problem  is  known  as  the  exact  solution.  The  difference  between  the  exact 
solution  and  an  approximate  or  hydrocode  solution  is  referred  to  as  the 
error.  For  purposes  of  comparison,  the  calculated  error  norms  are  the 
same  as  those  used  by  Hicks  and  are  defined  below. 

For  example,  take  a  variable  V  whose  exact  and  approximate  values 
in  cell  j  are  V  (j)  and  V  (j)  respectively.  Let  V  be  the  absolute  maximum 

of  V  (j)  where  j  ranges  over  the  region  of  interest.  For  V,  make  the 

E 

following  definitions: 

Sum.  Abs.  Error  =  h~  'VA^  "  VE^ 

M  ° 

Sum.  Sqr.  Error  =  -~  ?  (Vd)  "  Vj)) 

V    J 

M 

Maximum  Error  =  max  I  VA(J)  -  VE(J)1    sign  (V.(j  )  -  V  (j  )) 

V 
M 
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where  J,,  is  the  index  of   .   '   A       Ev 


For  each  test  problem  and  each  of  the  three  schemes,  error  tables  are 
given  containing  sum  abs.  error,  sum  sqr.  error  and  maximum  error.   Large 
values  of  the  maximum  error  norm  can  occur  at  the  shock  front  due  to  the 
numerical  approximations  of  the  shock  being  spread  over  two  or  three  cells. 
Thus,  the  maximum  error  norm  is  an  unreliable  norm  with  which  to  compare 
these  methods. 

For  all  test  problems,  a  polytropic  equation  of  state  is  used 
as  in  equation  (27). 

P  =   (7  -  1)  P  I  (27) 

For  all  test  problems,  7=  l.k. 

5.2  Test  Problem  1 
(a)  Description 

In  this  problem,  a  piston  moves  into  a  fluid  at  constant  velocity 
from  left  to  right.  The  fluid  is  consequently  compressed.  At  the  start  of 
the  problem,  a  shock  moving  from  left  to  right  exists  50  meters  ahead  of 
the  piston.  The  exact  solution  is  two  constant  states  separated  by  the 
shock  discontinuity  as  shown  in  Figure  6. 


3^ 


vT 


vr 


Figure  6.      Test  Problem  1 


V 


V 


V 


V 


X 


V 


sound  speed  to  the  left  of  the  shock 
sound  speed  to  the  right  of  the  shock 
pressure  to  the  left  of  the  shock 
pressure  to  the  right  of  the  shock 
density  to  the  left  of  the  shock 
density  to  the  right  of  the  shock 
velocity  to  the  left  of  the  shock 
velocity  to  the  right  of  the  shock 
piston  velocity 
shock  velosity 
initial  position  of  shock 


MC  where  M  is  a  positive  constant 

r  c 


X       =     initial  position  of  piston 
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(b)  Numerical  Values  for  1A 
M   =  1 

Ax  =  1  meter 

h  /  2 

P   =  10  dynes/ cm 

(■  .3 

p   =10   gm/cm 

V  =  0 
r 

V,  a:  1.18  x  105  cm/sec 

V  ~  2.08  x  105  cm/sec 

s  ' 

k  2 

P»  ~  3.^7  x  10  dynes/cm 

Q  o     -    2.30  x  10"  gm/cm^ 

X   =  51  meters  in  this  report  and  50  meters  in  [3] 

Right  boundary  at  300  meters 

X   =1  meter  in  this  report  and  0  meters  in  [3] 
P 

(c)  Boundary  Conditions 

On  the  right,  a  transmittive  "boundary  and  on  the  left,  a  piston 
starting  at  1  meter  and  moving  to  the  right  with  velocity  V  . 

(d)  Plot  Times   0.  and  .1  sec. 

(e)  Numerical  Values  for  IB 
The  same  as  1A  except 

M   =  100 

V»  ~  1.18  x  10'  cm/sec 

8        2 
P»  ~  1.68  x  10  dynes/cm 

C»  ~  6.26  x  10  cm/sec 

V  ~  1.42  x  10^  cm/sec 

s  ' 

Pa  2:  5.99  x  10   gm/sec 
(f)  Plot  Times  0.  and  10"  3  sec. 
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Computer  Time  =  100  Sec. 

Sum  Abs. 
Error 


Table  1.   Errors  for  Test  Problem  1A 
Problem  Time  =  10  J  Sec. 
Continuous  Rezone 

Iteration  Number  =  ^6k 


Pressure 
Density 
Energy 
Velocity 


.885 
.655 
.1+65 
.737 


Computer  Time  =  105  Sec. 

Sum  Abs. 
Error 


Pressure 
Density 
Energy 
Velocity 


.835 

.638 

.1+87 
.709 


Computer  Time  =   105  Sec. 

Sum  Abs. 
Error 


Pressure 
Density 
Energy 
Velocity 


1.031 
.886 
.658 

I.0I+5 


OIL 


Sum  Sqr.  Maximum  Position  of  Maximum 

Error  Error          Error 

.101  .285  current  shock  position 

.056  .201  current  shock  position 

.019  .122  current  shock  position 

.1^3  .326  current  shock  position 

Iteration  Number  =  36^ 

Sum  Sqr.  Maximum  Position  of  Maximum 

Error  '  Error          Error 

.076  .2^0  current  shock  position 

.OUl  .15*+  current  shock  position 

.018  .109  current  shock  position 

.133  .285  current  shock  position 

Iteration  Number  =  36^ 

Sum  Sqr.  Maximum  Position  of  Maximum 

Error  Error          Error 

.  LOU  -.236  current  shock  position 

.078  -.238  current  shock  position 

.027  -.115  current  shock  position 

.236  -.390  current  shock  position 


Donor 
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Computer  Time  =  99  Sec, 

Sum  Abs. 
Error 


Table  2.   Errors  for  Test  Problem  IB 
Problem  Time  =  10  J  Sec. 
Continuous  Rezone 

Iteration  Number  =  300 


Pressure 
Density 
Energy 
Velocity 


.9^5 
.963 
.980 
.828 


Computer  Time  =  98  Sec, 

Sum  Abs. 
Error 


Pressure 
Density 
Energy 
Velocity 


1.113 
I.165 
1.229 
1.005 


Computer  Time  =  98  Sec, 


Sum  Abs. 

Error 

Pressure 

1.569 

Density 

1.752 

Energy 

1.715 

Velocity 

1.370 

Sum  Sqr. 
Error 

.119 

.100 

.358 

.326 


OIL 


Sum  Sqr.   Maximum 
Error      Error 


Maximum     Position  of  Maximum 
Error  Error 

-.25^  current  shock  position 

-.266  current  shock  position 

-.587  current  shock  position 

-.556  current  shock  position 

Iteration  Number  =  299 

Position  of  Maximum 
Error 

.171  -.356  current  shock  position 

.157  -.366  current  shock  position 

.U63  -.658  current  shock  position 

.U^3  -.6Ul  current  shock  position 

Iteration  Number  =  300 

Position  of  Maximum 
Error 

.312  -.^95  current  shock  position 

.320  -.50^  current  shock  position 

.658  -.720  current  shock  position 

.6537  -.719  current  shock  position 


Donor 


Sum  Sqr.   Maximum 
Error      Error 
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(g)  Discussion  of  1A 

For  each  of  the  three  schemes  run,  there  are  spikes  at  the  shock 
front  and  slight  waves  behind  the  shock  front  as  shown  in  Figure  7.  Although 
the  waves  travelling  back  from  the  shock  front  are  larger  in  the  donor  method, 
the  donor  spike  obtained  at  the  shock  front  is  smaller  than  the  spike  for  OIL 
or  continuous  rezone. 

PUFF  has  very  good  shock  front  definition  and  slight  rarefaction 
dips  behind  the  shock.   Lax  Wendroff  and  PUFF  both  have  spikes  in  density 
and  internal  energy  at  the  initial  position  relative  to  the  piston,  i.e. 
50  meters  in  front  of  the  piston.   In  addition,  Lax  Wendroff  also  has 
spikes  at  the  shock  front  for  all  parameters. 


J\ 


Figure  7*   Typical  Wareform  for  Problem  1 
(h)  Discussion  of  IB 

Continuous  rezone,  OIL  and  donor  show  the  same  basic  behavior 
as  in  1A.  Again,  waves  travel  back  from  the  shock  front  as  in  1A.  The 
largest  waves  are  also  obtained  in  the  donor  scheme.   However,  the  small 
waves  may  be  preferable  to  the  very  large  spikes  obtained  in  PUFF  and  Lax 


39 


Wendroff .   PUFF  has  better  shock  front  definition  than  any  of  the  three 
schemes  tried  here  but  at  the  expense  of  large  spikes  50  meters  in  front 
of  the  piston.  Lax  Wendroff  has  oscillations  both  at  the  piston  face 
and  at  the  shock  front. 

5.3  Test  Problem  2 
(a)   Description 

In  this  problem,  a  piston  moves  with  constant  velocity  from 
right  to  left  thus  evacuating  a  pipe  which  initially  contains  constant 
state  gas  at  rest  as  in  Figure  8.   The  motion  of  the  piston  creates  a 
rarefaction  wave  moving  to  the  right. 


P    p 
r    r 


Figure  8.   Test  Problem  2 


(b)  Numerical  Values  for  2A 

h  /  2 

P   =  10  dynes/ cm 


-6     ? 
10   gm/cm 


C  =  1.4  x  10   cm  /sec 

r  ' 

V  =   -  C  /(y   +  1) 

p       r  ' 

Ax  =  100  cm 


Uo 


X  (o)   =  1300  meters  in  this  report  and  100  meters  in  [3]  . 

X      =  1500  meters  in  this  report  and  300  meters  in  [3]  . 

H 

(c)  Boundary  Conditions 

On  the  right,  a  transmittive  boundary  and  on  the  left,  a  piston 
moving  to  the  left  away  from  the  gas  with  velocity  V  . 

(d)  Plot  Times   .1  seconds 

(e)  Numerical  Values  for  2B 
the  same  as  2A  except 

V   =     r 


7+1 


Computer  Time  =  92  Sec. 


Table  3.   Errors  for  Test  Problem  2A 
Problem  Time  =  .1  Sec. 
Continuous  Rezone 

Iteration  Number  =  296  Sec. 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Posi 

tion 

of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.071 

.016 

-.031 

XK 

Density 

.821 

.011 

-.026 

XR 

Energy 

.456 

.003 

-.014 

XE 

Velocity 

2.1+7 

.111 

-.088 

x„ 

Computer  Time  =  93  Sec. 

Sum  Abs. 
Error 


Pressure 


1.096 


OIL 


Iteration  Number  =  296 

Sum  Sqr.    Maximum    Position  of  Maximum 
Error       Error  Error 

.017     -.032  X_ 


4l 


Density- 

.848 

.011 

-.026 

Energy 

.470 

.003 

-.014 

Velocity 

2.536 

.115 

Donor 

-.Q89 

Computer  Time  =  90  Sec, 


LE 


XT 


X 


R 


Iteration  Number  =  296 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

3.102 

.218 

.135 

rarefaction  wave 

Density 

2.726 

.151 

.120 

rarefaction  wave 

Energy 

1.088 

.019 

.042 

rarefaction  wave 

Velocity 

5.586 

.680 

.269 

rarefaction  wave 

Table  4.   Errors  for  Test  Problem  2B 
Problem  Time  =  .1  Sec. 
Continuous  Rezone 
Computer  Time  =  95  Sec. 


Iteration  Number  =  296 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.226 

.014 

-.023 

2600  metres 

Density 

1.052 

.011 

-.023 

2600  metres 

Energy 

.924 

.007 

-.017 

2600  metres 

Velocity 

2.034 

.055 

-.057 

2600  metres 

OIL 


Computer  Time  =  101  Sec. 

Sum  Abs. 
Error 


Sum  Sqr. 
Error 


Maximum 
Error 


Iteration  Number  =  296 

Position  of  Maximum 
Error 


Pressure 


1.201 


.013 


.023 


2600  metres 


1+2 


Density 

1.021+ 

.011 

-.023 

2600  metres 

Energy 

.906 

.006 

-.017 

2600  metres 

Velocity 

1.998 

.053 

Donor 

-.055 

2600  metres 

Computer  Time  =  96  Sec. 


Iteration  Number  =  296 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

2.818 

.172 

.135 

2786  metres 

Density 

2.528 

.120 

.122 

2786  metres 

Energy 

1.101+ 

.017 

.01+2 

2788  metres 

Velocity 

3.072 

.164 

.133 

2778  metres 

(f )  Discussion  of  2A 

The  results  obtained  for  continuous  rezone  and  OIL  are  very 
good.   Donor  displays  the  start  of  strong  instabilities  in  the  rarefaction 
wave.   In  continuous  rezone,  OIL  and  donor,  there  is  a  slight  rise  "in 
density  and  slight  drop  in  specific  internal  energy  just  in  front  of  the 
piston  as  in  Figure  9  which  seems  to  have  little  or  no  effect  on  the  rare- 
faction region.  This  effect  is  slightly  worse  in  donor  than  the  other 
two.   The  reason  for  this  effect  might  be  connected  with  always  choosing 
the  donor  value  of  density  in  the  density  weighting.  Also,  OIL,  continuous 

rezone  and  donor  underround  at  X_  and  at  X  .   PUFF  and  Lax  Wendroff  show 

L  K 

undershoots  to  the  left  of  X  .   Lax  Wendroff  even  oscillates  slightly  at 

K 

the  left  of  X_. 
K 


**3 


Xr         Xc 

Figure  9-   Internal  Energy  for  OIL  and  Continuous  Rezone 

(g)  Discussion  of  2B 

The  "basic  behavior  of  2B  is  similar  to  that  of  2A.  Again,  donor 

breaks  down  in  the  rarefaction  wave.  PUFF  and  Lax  Wendroff ,  due  to  the 

undershoots  obtained  to  the  left  of  X^,  seem  inferior  to  continuous  rezone 

K 

and  OIL  for  this  problem. 


5.U  Test  Problem  5 
(a)  Description 

This  problem  is  a  special  case  of  the  Riemann  problem  and  is 
known  as  either  the  shock  tube  problem  or  the  ruptured  membrane  problem. 
The  problem  is  to  find  the  hydrodynamic  variables  after  the  removal  of  a 
membrane  at  time  t  =  0  separating  two  constant  unequal  states  at  rest  as 
in  Figure  10.  We  adopt  the  convention  P  >  P  >  0  and  investigate  the 

three  cases  pp>p  >p0    =  p  ,p,<p  .   When  the  membrane  has  been  removed, 

*   *  *    r  Hi       r 

a  rarefaction  wave  travels  into  the  left  state  and  a  shock  travels  into 
the  right  state. 


pi    pi 


hk 


membrane 


P     p 
r     r 


X 


Figure  10.   Test  Problem  5 


(b)  Numerical  Values  for  5A 

X  (0)   =  101  metres  in  this  report  and  100  metres  in  [3] 


Ax 

P 


I 


I 


V 


V 


=  1  meter 

8        2 
=  10  dynes/cm 

=  10"5  gm/cm3 

=  0 

h  ,  2 

=  10  dynes/cm 

=  10   gm/cm 


=  0 


Left  boundary  at  2  metres  in  this  report  and  0  metres  in  [3]  . 
Right  boundary  at  252  metres  in  this  report  and  250  metres  in  [3]  • 

(c)  Boundary  Conditions 

Transmittive  boundaries  exist  at  both  the  left  hand  and  right 
hand  boundaries. 

(d)  Plot  Times  2  x  10  J  sec. 


h<? 


(e)  Numerical  Values  for  5B 

The  same  as  5A  except 

X  (0)  =  251  meters  in  this  report  and  250  meters  in  [3] 
s 

(■  Q 

p      =10   gm/cm 

Right  boundary  at  502  metres  in  this  report  and  500  meters  in  [3]« 

(f )  Numerical  Values  for  5C 

The  same  as  5A  except 
X  (0)  =  251  meters  in  this  report  and  250  meters  in  [3] 

p      =10   gm/cm 

-5      3 
p      =10   gm/cm 

Right  boundary  at  502  metres  in  this  report  and  500  meters  in  [3]. 


Table  5»   Errors  for  Test  Problem  5A 
Problem  Time  =  .002  Sec. 
Continuous  Re zone 


Computer 

Time  =  85 

Sec. 

Iterat 

ion 

Number  =  2< 

Sum  Ab£ 

<  • 

Sum  Sqr. 

Maximum 

Posit 

ion 

of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.480 

•037 

-.123 

xs 

Density 

3.199 

•293 

-.3^8 

xs 

Energy 

2.890 

.2kk 

-.273 

xs 

Velocity 

3-357 

.927 

-.832 

x= 

1+6 


OIL 


Computer  Time  =  86  Sec 


Iteration  Number  =  200 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position 

of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.1+81 

.01+1 

-.137 

xs 

Density- 

3.292 

.31+1+ 

-.392 

xs 

Energy 

2.911 

.251+ 

-.278 

xs 

Velocity 

3- 518 

1.088 

Donor 

-.861 

xs 

Computer 

Time  =  88  Sec. 

Iteration  Number  =  ; 

Sum  Abs. 

Sum  Sqr. 

Maximum 

Position 

of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.771 

.072 

-.11+0 

xs 

Density 

3.602 

.1+15 

-.1+03 

xs 

Energy 

2.801+ 

.268 

-.277 

xs 

Velocity 

3.630 

1.229 

-.861 

Xo 

(g)  Discussion  of  5A 

This  problem  demonstrates  very  clearly  the  instabilities 
which  can  occur  in  rarefaction  regions  with  the  donor  scheme.   The 
instabilities  in  this  particular  problem  seem  to  be  caused  by  the 
slope  to  the  right  of  X  on  all  donor  plots  being  too  steep  as  shown 
in  Figure  11.  A  general  loss  of  fine  detail  is  also  in  evidence  in 
all  three  methods.   In  particular,  there  is  strong  underrounding  of 
the  density  at  X^  and  a  consequent  loss  of  definition  at  X^  in  the 
specific  internal  energy  plot.  The  continuous  rezone  and  OIL  plots 
are  essentially  identical  but  for  minuscule  peaks  at  X  on  the  pressure 


hi 


and  velocity  in  continuous  rezone.   PUFF  causes  less  smearing  of  the 
density  at  X^  than  in  OIL  or  continuous  rezone.   However,  Puff  also 

has  overrounds  at  X  which  are  not  in  evidence  with  continuous  rezone 

K 

and  OIL.   The  Lax-Wendroff  solution  has  several  spikes. 


Table  6.  Errors  for  Test  Problem  5B 
Problem  Time  =  .002  Sec. 
Continuous  Rezone 
Computer  Time  =  172  Sec. 


Iteration  Number  =592 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Po 

sition  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1-358 

.062 

-.218 

1+01  meters 

Density 

7-758 

2.0&U+ 

•  ^599 

376  meters 

Energy 

9-588 

H.136 

.61+9 

375  meters 

Velocity 

3-137 

•  679 

-.736 

1+01  meters 

OIL 


Computer  Time  =  176 


Iteration  Number  =  592 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.391 

.072 

-.239 

1+01  meters 

Density 

7-882 

2.157 

-.1+87 

1+01  meters 

Energy 

9-5^9 

1+.103 

.61+8 

375  meters 

Velocity 

3.2U1 

•  753 

Donor 

-.760 

1+01  meters 

Donor  will  not  run  to  .002  sec.  due  to  instabilities. 


hQ 


Xo     Xs 


Figure  11.   Donor  Density  in  Problem  5A 

(h)  Discussion  of  5B 

The  instabilities  in  the  rarefaction  region  for  the  donor 
scheme  are  so  bad  in  5B  that  the  problem  will  not  run  to  .002  seconds. 
Even  at  .001  seconds,  the  instabilities  in  donor  are  quite  bad.   As  in 
5A,  continuous  rezone  and  OIL  are  essentially  identical.   Also,  the 
general  behavior  of  OIL  and  continuous  rezone  is  similar  to  that  in 
5A.   In  PUFF,  the  overshoot  in  velocity  is  quite  pronounced,  although 
the  definition  at  X_  is  still  good.   The  Lax  Wendroff  solution  shows 


oscillations. 


h9 


Table  7.   Errors  for  Test  Problem  5C 
Problem  Time  =  .002  Sec. 
Continuous  Re zone 


Computer  Time  =  178  Sec. 


Iteration  Number  =592 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1-3^3 

.11+0 

-.328 

312  meters 

Density 

6.356 

2.206 

•571 

302  meters 

Energy 

12.111 

8.737 

.890 

301  meters 

Velocity 

i+.oi+o 

.655 

OIL 

-.699 

312  meters 

Computer  Time  =  176  Sec. 


Iteration  Number  =592 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.39*+ 

.11+8 

-.31+0 

312  meters 

Density 

6.397 

2.233 

•  575 

302  meters 

Energy 

12.101 

8.717 

.890 

301  meters 

Velocity 

4.208 

.682 
Donor 

-.705 

312  meters 

Donor  will  not  run  to  .002  sec  due  to  instabilities, 
(i)  Discussion  of  5C 

As  in  5B,  donor  will  not  run  to  .002  seconds  due  to  severe 
instabilities  in  the  rarefaction  region.   Overrounding  occurs  with 
both  continuous  rezone  and  OIL  in  the  pressure  and  velocity  at  X^. 
The  overrounding  seems  slightly  worse  in  OIL.  However,  the  overrounding 
at  X^  in  the  PUFF  solution  seems  worse  than  either  continuous  rezone 
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or  OIL.   The  handling  of  the  density  at  X^  by  PUFF  is  considerably 
better  than  continuous  rezone  or  OIL.   The  Lax-Wendroff  solution  is 
very  poor  with  many  overshoots  and  spikes  etc. 

5-5  Test  Problem  6 
(a)  Description 

This  problem  is  another  special  case  of  the  Riemann  problem 
and  involves  the  collision  of  two  shocks  as  in  Figure  12.  The  first 
part  of  the  problem  investigates  the  case  when  P«>P  •  The  second  part 
examines  the  collision  of  two  symmetrical  shocks. 


pl 

p 

r 

pr 

1 

■ ■■           

V 
r 

X 


s^(t) 


X 


Sr(t) 


Figure  12.   Test  Problem  6 

(b)  Numerical  Values  for  6A 

Ax  =  1  meter 

X  -(0)  =  75  meters 

X  (0)   =  125  meters 
srv  ' 

P   =  10  dynes /cm 
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p   =10   gm/cm 


V 


0  cm/s 


ec 


Q  o 

P   =  10  dynes  /cm 

7         2 
P   =10  dynes/cm 

r  ' 

Left  boundary  at  70  metres  in  this  report  and  0  metres  in  [3].   Right 
boundary  at  180  metres  in  this  report  and  200  metres  in  [3]. 

(c)  Boundary  Conditions 

Transmittive  boundary  conditions  are  applied  at  both  the 

left  hand  and  right  hand  boundaries. 

-4 

Plot  Times  0.,  and  7  x  10   seconds 

(d)  Numerical  Values  for  6B 

The  same  as  6A  except 

O  p 

P  =  P   =10  dynes/cm 

r  ' 


Table  8.   Errors  for  Test  Problem  6A 

Problem  Time  =  .0007  Sec. 

Continuous  Rezone 


Computer 

Time  =  75 

Sec 

Iteration  Number  =  . 

Sum  Ab£ 
Error 

• 

Sum  Sqr. 
Error 

Maximum 
Error 

Posit 

ion  of  Maximum 
Error 

Pressure 

1.477 

.431 

•  515 

-* 
XSR 

Density 

1.914 

•  373 

.412 

-* 
XSR 

Energy 

1.924 

•  334 

.336 

x  1 
col 

Velocity 

1.104 

.148 

-.229 

XS1 
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OIL 


Computer  Time  =  75  Sec, 


Iteration  Number  =  163 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position 

of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.298 

.309 

•398 

XSR 

Density 

1.784 

•293 

•307 

-* 
XSR 

Energy 

1.866 

•  314 

•  330 

x  1 
col 

Velocity 

1.064 

.120 

Donor 

-.226 

* 

XS1 

Computer 

Time  =  79  Sec. 

Iteration  Number  =  ! 

Sum  Abs. 

Sum  Sqr. 

Maximum 

Position 

of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1.315 

.369 

.465 

* 
XS1 

Density 

2.166 

.449 

-.397 

* 

XSR 

Energy 

1.951 

.325 

•297 

x  1 
col 

Velocity 

1.375 

.305 

-.443 

XQT3 

(f )  Discussion  of  6A 

Continuous  rezone  and  OIL  both  have  small  overshoots  at 
X„  in  pressure,  density  and  velocity  as  shown  in  Figure  13  while  donor 
has  almost  no  overshoot.  Also,  there  is  a  general  smearing  of  fine 
detail  as  has  been  noticed  in  previous  problems.  The  specific  internal 
energy  plots  for  continuous  rezone  and  OIL  are  slightly  better  than 
donor.   OIL  and  continuous  rezone  are  practically  identical.  PUFF 
has  large  spikes  in  the  density  and  specific  internal  energy.  The 
Lax-Wendroff  solution  is  barely  recognizable  due  to  large  oscillatory 
components. 
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XsR 


Figure  13-  Pressure  for  Problem  6A 


Table  9.  Errors  for  Test  Problem  6B 
Problem  Time  =  .0007  Sec. 
Continuous  Rezone 


Computer 

Time  =  73  Sec. 

Iteration  Number  =  . 

Sum  Abs. 

Sum  Sqr. 

Maximum 

Po£ 

sition  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

5.51+1 

1.336 

.1+93 

83  metres 

Density 

1+.322 

.808 

.31+7 

83  metres 

Energy 

2.571 

.1+29 

.382 

83  metres 

Velocity 

30.831 

26.560 

-.990 

76  metres 

51+ 


OIL 


Computer  Time  =  75  Sec 


Iteration  Number  =  162 


Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

1+-923 

1.163 

.1+63 

83  metres 

Density 

3.91+8 

.69I+ 

.312 

83  metres 

Energy 

2.308 

.381+ 

.368 

83  metres 

Velocity 

31.001 

26.I+9I+ 

Donor 

-.990. 

76  metres 

Donor  -will  not  run  to  .0007  sec.  due  to  instabilities. 


-J 


Figure  ll+.   Pressure  for  Continuous  Rezone  in  Problem  6B 
(g)  Discussion  of  6B 

More  difficulty  was  obtained  with  this  problem  than  any 
other.  The  donor  scheme  for  6B  would  not  run  to  7  X  10~  seconds. 
The  reasons  for  this  difficulty  are  not,  as  yet,  understood.  It  can- 
not be  blamed  on  rarefaction  regions,  since  none  exist  in  this  problem. 
The  donor  scheme  runs  up  to  and  just  beyond  the  collision  but  then 
'blows  up'.   Both  OIL  and  continuous  rezone  show  serious  oscillations 


55 


soon  after  collision  as  in  Figure  1^.   However,  no  asymmetries  were 
observed.   Several  variations  of  the  basic  continuous  rezone  scheme 
were  tried  not  using  artificial  viscosity,  but  none  of  them  showed 
fewer  oscillations  than  the  standard  method.   PUFF  shows  some  large 
spikes  in  the  density  and  specific  internal  energy  but  the  pressure 
and  velocity  are  very  good.   The  Lax  Wendroff  solution  has  instabilities 
comparable  to  OIL  and  continuous  rezone.   In  summary,  the  collision  of 
two  asymmetric  shocks  works  reasonably  well  but  a  lot  of  difficulty  is 
obtained  with  the  collision  of  two  symmetrical  shocks.   From  the  plots, 
it  can  be  seen  that  the  instabilities  for  continuous  rezone  and  OIL 
arise  in  the  region  where  the  velocity  is  supposed  to  be  zero.   Hence, 
in  this  region  the  effective  viscosity  will  be  close  to  zero  and  there 
will  be  almost  no  damping  effect  for  any  slight  perturbations. 

5.6  Test  Problem  7 
(a)  Description 

This  problem  is  concerned  with  the  overtaking  of  one  shock 
wave  by  another  as  in  Figure  15-  When  two  shock  waves  are  travelling 
in  the  same  direction,  the  one  behind  will  always  overtake  the  one  in 
front.   This  is  a  consequence  of  the  fact  that  the  front  side  of  a 
shock  is  supersonic  relative  to  the  shock  front  and  the  backside  is 
subsonic.   When  the  overtaking  has  occurred,  a  rarefaction  wave  travels 
to  the  left  and  a  stronger  shock  continues  moving  to  the  right. 
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Figure  15-   Test  Problem  7 

(b)  Numerical  Values  for  7  A 

Ax  =  1  meter 

h  ,  2 

P   =  10  dynes/ cm 

p   =10   gm/cm 

V   =  0 
r 

-,^8      /  2 

P.    =  10  dynes/ cm 

12         2 
P   =  10   dynes/cm 

X   (0)   =  25  meters 

X   (0)   =  100  meters 
sr 

Left  boundary  at  2  meters  in  this  report  and  0  meters  in  [3].  Right  boundary 
at  200  meters . 

(c)  Boundary  Conditions 

Transmittive  boundary  conditions  are  applied  at  both  the  left 

hand  and  right  hand  boundaries. 

-5 

(d)  Plot  Times  0.  and  3  x  10   seconds 
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Table  10.  Errors  for  Test  Problem 
Problem  Time  =  3.  X  10  J   Sec. 
Continuous  Rezone 


7A 


Computer 

Time  =  98  Sec. 

Ite 

ration  Number  = 

Sum  Abs. 

Sum  Sqr. 

Maximum 

Pos 

ition  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

2.667 

.167 

.142 

185  meters 

Density 

3.230 

.208 

-.135 

172  meters 

Energy 

5.146 

1.229 

.U65 

172  meters 

Velocity 

1.581 

.181 

-.340 

186  meters 

OIL 


Computer 

Time  =  99  Sec. 

Iteration  Number  =  ' 

Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

2.772 

.169 

.139 

185  meters 

Density 

3.260 

.207 

-.138 

172  meters 

Energy 

5.376 

1.316 

.474 

172  meters 

Velocity 

I.687 

.214 

Donor 

-.374 

186  meters 

Computer 

Time  =  93  Sec. 

Iteration  Numbe] 

Sum  Abs. 

Sum  Sqr. 

Maximum 

Position  of  Maximum 

Error 

Error 

Error 

Error 

Pressure 

3.024 

.177 

.155 

185  meters 

Density 

3.330 

.202 

-.140 

172  meters 

Energy 

5.73^ 

1.430 

.485 

172  meters 

Velocity 

1.886 

.260 

-.397 

186  meters 

=  331 
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(e)  Discussion  of  7  A 

_5 
At  3  X  10   seconds,  there  is  a  slight  dip  to  the  left  of  the 

rarefaction  wave  on  all  plots  for  each  of  donor,  continuous  rezone  and 

OIL  at  ninety  meters  as  in  figure  l6.   One  does  not  know  if  these  dips 

are  obtained  with  PUFF  and  Lax  Wendroff  since  in  Hicks  report,  the  plots 

for  this  problem  start  at  115  meters.  A  general  smearing  of  shock  corners 

is  evidenced  as  in  other  problems.   Continuous  rezone  has  tiny  overshoots 

in  the  pressure  and  density  at  X  which  do  not  arise  in  OIL  or  donor.  The 

s 

plots  obtained  for  velocity  and  specific  internal  energy  for  continuous 
rezone,  OIL  and  donor,  are  practically  identical.  Also,  there  is  very 
little  difference  between  continuous  rezone  and  OIL  in  general.   PUFF 
has  good  overall  definition  but  there  are  two  very  large  spikes  on  the 
density  and  specific  internal  energy. 


90  METERS 


Xs 


Figure  16.  Pressure  for  Continuous  Rezone  in  7A 
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6.   GENERAL  DISCUSSION 

A  comparison  of  the  three  differencing  schemes  described 

[3] 

in  this  report  with  PUFF  and  Lax  Wendroff    seems  to  refute  the  often 

asserted  statement  that  Lagrangian  codes  are  overwhelmingly  superior 

[12] 
to  Eulerian  codes  in  one  dimension. 

In  general,  donor  preserves  quite  good  definition  at  the 
shock  front  but  is  very  unstable  in  rarefaction  regions.  Donor  also 
shows  waves  behind  the  shock  front  (on  the  rarefaction  side)  which  are 
evidence  of  the  start  of  instabilities.   Continuous  rezone  and  OIL  are 
stable  in  rarefaction  but  often  have  larger  overshoots  than  donor  at 
the  shock.   For  many  of  the  test  problems,  continuous  rezone  and  OIL 
are  almost  indistinguisable.   There  is  a  general  loss  of  fine  detail 
in  rarefaction  regions  in  all  three  schemes  as  is  evidenced  by  problems 
5  and  6A. 

The  mechanism  for  breakdown  of  donor  in  rarefaction  is 
not,  as  yet,  ascertained.  In  problem  5,  instabilities  seem  to  set 
in  at  the  intersection  of  two  curves  of  differing  gradients,  whereas 
in  OIL  and  continuous  rezone  there  is  one  continuous  curve. 

For  a  general  purpose  hydrocode,  donor,  due  to  its  rarefac- 
tion problems  is  less  preferable  than  either  OIL  or  continuous  rezone. 
A  hybrid  scheme  which  uses  donor  at  the  shock  front  and  one  of  the 
other  two  codes  in  the  rest  of  the  problem  is  probably  best. 
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APPENDIX  A 
PLOTS  FROM  THE  TEST  PROBLEMS 
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APPENDIX  B 
DERIVATIONS  OF  ORDERS  OF  TRUNCATION  ERRORS 

Conservation  of  Mass 

(a)  Continuous  Rezone 

It  is  necessary  to  find  expressions  for  U.  1 

X~2 

From  equation  (12) 


un  =  un  -  (p.nn  -  p.nJ  -&. 

1      i     i-fl    i-ly  _  n 


Therefore  Un  =  U? 
l      l 


2p.  Ax 

l 


fx   |  ,  .  o(ax2)]  &  •    (IT) 

pi 


Hence  U.nl  =  U.nl  -  —  |£  I  .  +  O(AxAt)  (l8) 

i 

noting  that  r-=  ~  1  +  E    if  E  «1 
JL  —ill 


Also  t)A  =  UA  .  *  ||  |   +  O(AxAt)  (19) 

Therefore 

!J.*1  -  ffA  =  g-  |  .  ax  +  0(ax3)  +  O(AXAt)  (20) 


n       n   *n 
In  general  AM,  1  =  p.  1  U.  1  At 

2      1+2  1+2 


Ik6 


Hence   AM.  1 
1+2 


11  TT   nT    A  4- 

p.  tf1+l  At 


1  +  (u.ni  -u.ni)  ^ 

1+2    i~2  Ax 


pn  U  A  At 
1   1+2 


1  +  |2  I  .  At  +  0(At2)  +  0(Ax2At) 

OX  '  1 


from  equation  (20) 


Therefore 

A,,  nn      n  ., 
AM .  ,  1  =  p  .  At 

1+2 


U  A  -  At 
1+2 


_1  dP 
n  Sx" 


au 


.  +  U* 

i    i  ox 


+  0  (AxAt)  +  0(At  ) 


(21) 


from  equation  (19) 


Also  AM.nl 
X"2 


11      AX 

Pi-1  ^ 


U.ni  -  At 
1"2 


n  ox  '  i    l-l  ox 


i-1 


+  O(AxAt)  +  0(At  ) 


The  equation  of  continuity  is 


Mn+1  -  M*  =  AM.ni   -  AM.nl 


i+: 


Therefore 


At 


i   =  P.  ,  U.  1  -  p.  U.  1  +  0(At  )  +  O(AxAt) 
Ki-1  1-2    i   1+2 


lU7 


=  P^!  U^!  "  P?  lI?  +  ^A*2)  +  °(At2)  +  O(AxAt) 


Therefore 


Pn+1  -  Pn    (Pu)n  -  (Pu)n  - 
1      1     __  1      1-1 


At 


Ax 


-  It +  ^  +  °C*0  ♦  °(a*)  +  <  ) 


(b)   OIL 


AMA 
1+2 


pnAt  U.ni 
1+2 


!  t  (u*   -Un)  ^ 
v  l+l    iy  Ax 


(22) 


pn  At  Unl 
1+2 


1  +  ^7  I  niAt  +  O(A^At)  +  0(At2) 
ox   1+?T 


Therefore 

AM.nl  =  pn  At 
i+2    x 


U.ni  -  At 
1+2 


1_  dP 

n  d~x 
P^ 


du 


l    ■  l  dx 


+  O(AxAt)  +  0(At  ) 


(23) 


Since  equation  (23)  is  identical  to  equation  (2l)?  the  rest  of  the  calculation 

for  OIL  is  the  same  as  for  continuous  rezone. 

(c)   Donor    AM.nl  =  p."  U."  At 
1-7T      i-l   i-l 


Also  p 


n  ~n     n  n     P    -  P.  A 
.  U.   =  p.  U   -   i+1    i-l 
l   l      l   l     L  -J 


At 
2Ax 


from  equation  (12) 


11+8 


Therefore 

P^1  -Pn  Pnun  -P»  U*         ,2 

i i  ,  1   1     1-1   l-l   _    ,  d   p   | 

aT"+  S         At^2  'i-i 


Therefore 

n+1    n     n  TTn     n  TT  n 
p.    -P.    p.  IT.   -  p.  .  U.  .    ^    >,  TT 


Conservation  of  Momentum 
(a)  donor 
The  form  of  the  conservation  of  momentum  equation  which  will  be  used  is 


dpU   dpU2   dP    _ 
3T+  37"  +  Sx"  =  ° 


The  finite  difference  equation  used  is 


M?+1  tf"1       =     M^U*     +  AM.nl     U.    ,    -AM.nl     U.  (2k) 


i-p        i-1  i+o        i 


Therefore 


M?+1U?+1      -     I^U"  '        „       „       „ 

l_i.     p."    Tj.     2.  pnu.2     -     (P."     -  P.M 
i-1     l-l         i     i  l+l         i-l 


At 


Therefore 

n+1  TTn+-l         n  TTn  _  n         _,  n               „         ~  2             ~  2 

P .        u .        -  p .    u.  P.  , ,    -  P.    -              p .    t   U.    ,    -  p  .U.                        /oc-\ 

i          j l     l  l+l         i-l     _         i-l     i-l         l  l                       (25) 

At  2  Ax              ~                         Ax 


li+9 


The  right  hand  side  equals 


pi-l 


l-l 


i 2  r  i 

f^P  I  n^A  2\\  £b  n     tt11      Ap  I  ft^2^At 


Ax 


from  equation   (17) 


n      /._  n   s2         n   /T7n\2 
pi-l   (Ui-l>      -  pi    (U1} 

Ax 


+   O(At) 


Therefore 


>n+i  un+l  -  Pn  u2?  Pn  d^)2  -  P.nn  (uny  p.*  -  p.n 
1    1     ki  1    1 v  1      1-1 v  1-1      1+1    1- 


At 


Ax 


2Ax 


|pu  +  Ipu2  +  |p  +  Q(Ax)  +  Q(At) 


dt  dx  3x 


(t>)   continuous  rezone 


Un  AM.    1  =  pn  At 
1        1+2  1 


u.ni 


n  dx  I  1  1  ox  I  1 

pi 


_+  O(AxAt)   +  0(At    ) 


dP  .  rt/A  2v\  At 

^ii  +  0(Ax  n- 


pi 


1 


from  equations    (17)   and  (2l) 
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i         i 


At 


n  AM.ni  ff.nn    -  AM.    1  9? 

,  1-7T      1-1  1+7T      1 

i i  i+l         i-1  2  2 


pn+l  ^+1     _  pn  yn       p_^   _  p^ 


2Ax 


AxAt 


from  equation   (25) 


The  right  hand  side  equals 


'i-1 
Ax 


u.ni  u.nn 

i-^  i-i 


-  At  U 


n   I  -i- 
H'i-i 


1       dP  ,  „n     SU 


ox  '  1  i-i  ox  '  1 


-'} 


-  Ui^(ii  I  i  i +  0(Ax2))  ^r +  °(^t)  +  0(At2) 

2  v.  J  p  J 


n 
Ax 


U.nl  U"  -  Atl^ri-^l.    +  U?|H 
i+7T     i  i  (     n  ox  '  l  l  ox 

lpi 


) 


"  Ui^ff  I  i  +  °(^2))^  +  O(AxAt)   +  0(At2) 


1    j      n     TT  nn  TT  n  n  TT  n.   TTn\ 

—     p.    t   U.    1  U.    .    -  p .   U.    1  u.  ) 
Ax  \    i-l     i~2     1"1         x     1+2     V 


r    TT  n  j  oP  ,  n     TT  n     oU  |        \ 

1-13  ox  '  i  i-l     i-l  ox  '  i-l/ 


+  0(At)   +   0(||   ) 


At 
Ax 


1  \  OX      1  1      1   ox 


'} 


;«4(lii..-»'-2i}-»4(£'.*°to8| 
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^  +   O(Ak)   +   O(At) 


vAx 


Hence,    equation    (25)   for  continuous   rezone  becomes 


n+1  TTn+l         n  n  n  /TTns2            n  ,  _  n   ,2  _  n         _  n 

p.        U.        -  p.  IT.  p.  (U.)      -  p.    .  (U.    , )  P.  ,_    -  P.    . 

1           i 1  1  1  i              'l-l  v    i-l/  14-1           1-1 

At  Ax  +    2"Ax 


dpU  ,  SpU2,  dP  ,  ^/   v    _/Ajv    ^/At2 


at 


-  ^  +  §£■  4-  0(ax)  +  O(At)  +  0(^  ) 
dx    ox  VAX 


(c)   OIL 

Since  equation  (23)  is  identical  to  equation  (21)  and  the  U. 
(i  an  integer)  for  continuous  rezone  and  OIL  are  equal,  the  calculation 
for  OIL  is  the  same  as  for  continuous  rezone. 

Conservation  of  Energy 

The  form  of  the  conservation  of  energy  equation  which  will  be 

used  is 

SfiE  +  SpJU  +  aPU  =   Q 
~dt         dx  o\x 

It  can  be  seen  that 


152 


En 

l 


=   E 


p.  At 

1 

P.  Ax  h 
UnAt 


[u.n,  +u.nn  -  u.n,  -  u.nJ 

i+j.    l+l    1-1    i-lJ 


_  u  At  ,  n     n  s  d  ,     n 

— ^i+1    i-l/    At  ^  ^  lii+l 


2p*Ax 


8(p^)2Ax2 


^)2 


The  conservation  of  energy  equation  in  finite  difference  form  is 


T,/rn+l     -nn+l  r/P-     T=in  a*,     nn  ^     n  *»,     ^,      ~n 

M.        E.        =  K.    E.      +     AM.    1     E.    ,    -     AM. IE. 
li  li 


i~2        i"1 


1+2 


.,n+l  -n+1  -Ji  _n 

M.         E.         -  M.    E. 

or     _i l  ii 

At 


-F 


n 


rTT  n         ~  n         T.  n         ~  n   -, 
i      [U. , ,    +  U.    ,    -  U.    ,    -  U.    ,  1 
-j7         i+1  l+l  i-l  i-l 


-Ui    (P  n   -  P  n   ) 
-|  V±i+1     H-l; 


n  N2 


..n-t-1  _n+l       „n     n 

M.        E.         -  M.    E. 

Hence       i  l l     l 

At 


+_^_    (p^  -  p^) 

8pn  Ax 


n       *n,  ~  n  n     *n.,     ^n 

+  p .    _   U.    1  E.    ,    -  p.    U.  ,1     E. 

Mi-1     i~2     i-l       Mi      1+2        x 


-dPU    ,n   .      J        n       *n.,  ~  n 
•s —     .   Ax  +  p.    .    U.    IE.    , 
OX       '  1  'l-l      i-^-      1-1 


n     *nn     <^n 

-  p .    U.    1     E. 

Mi     1+2       i 


+   0(ax    )   +   O(AxAt) 


(26) 


(a)     Donor 


U*nl     =     u"     and  U*nl     =  Un 
i"2  i-1  1+2  x 


n       *n,  ~  n  n     *n..     ~n 

Hence     p .    _    U.    1  E.    ,    -  P .    U.    1     E.    = 
l-l     i~2     i-l  i     1+2        i 
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i-l     i-l 


-,1 


»  n 

^  n  P."     At 

Ei-1  "     -i=i (U?  +  un  -  u.n0  -  u.nQ) 

p.^Ax^       X  X  X"2  X"2 


U.n,    At  _ 

-1=1 —  (Pn  -  p.*)  +■  o(At2) 
n   A      v    l  i-2  / 


2p        * 


n  TTn 
piUi 


T-.n  _n   A, 

E.    -     P.  At  /..  n  ~  n 

1             l  (U.    ,  +  U.    ., 
-  v    l+l  i+l 

pnAx   1| 


-u"    -u.nJ 

i-l  i-l 


U?  At  0 

"  — (P*      -   P.nJ    +   0(At2) 

_  n  A      v   i+l         i-l'  v        ' 

2p  .   Ax 
l 


=     p.n_      U.n,      E.n.       -  pn  Un  En  +   0(At2)    +   O(AxAt) 
i-l        i-l        i-l  ill 


Hence  equation   (26)  "becomes 


n+l  „n+i         n  _n  ,  _Tvn       ,  _,T\n 

P.        E.        -  p.    E.  «v_T  (pEU).    -    (pEU).    . 

li  li         oPU    in  4-      v        ' i        x        /i-. 


At 


dx      'i 


Ax 


=     O(At)    +   O(^) 


15k 


Therefore 
,n+l 


(pE)f L  -   (pE)J  (  pEU)*   -   (PEU)^  (PU)J  -   (PU)^ 


At 


A  x 


Ax 


dpE         dpEU       dPU 
dt  chc  dx 


*   O(At)    +   O(Ax)    f   0(^-) 


(t>)   continuous  rezone 


n     "*n     ~n 

p.    U.    1  E. 

i      1+2     1 


n 


1+2  \,2.  P  J31     i+2         1+2  5x-|i| 


E^   -  p?  At      ...  n         ~  n 

1       -i (Vl  +  Ui+1 

pnAx   U 
l 


UnAt 
l 


n        ~  n 


-U."     -  U."   ) 

l-l  l-l 


^l"Pi-l)      +  °^> 


2PiAx 


Hence,    the  right  hand  side  of  equation   (26)  becomes 


.  |2L   |J  ax  +   0(Ax2)    +   O(AxAt) 


U.n     +|U   |  Ax     _At 

l-l       ox    'i-l       p 


1       dP 


-  v-  I  •    1     +  U.    1 


nn   dU 


i-p-  §x    '  i-li 


+   O(AxAt)    +   0(At2)    +   0(Ax2) 


X 
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P  *     At 
p  n        n     _  _i-l         [xjn  +  ^n  _  ,n     _  g  n  -, 

l-l     l-l  4Ax       l  l  1-2  i-2 


U 

•    (Pn   -   P.no)   At  +  0(At2) 


2Ax    v    i  i-2' 


TTn       dU  I       Ax        . .   f  1     dP  I        '        TT  n.    ^U  p 

U .     +   -T-       .    -7T    -    At  /  T-       ,1+    U.     It-      4 

i       ox     i     2  <  2  n  dx     i+^         i+^  dx     i 


+   O(AxAt)   +  0(At2)   +  0(Ax2) 


X 


PnAt 
pnEn._i  ^n         ~n  n    .  jjn    j 

l      l  4Ax        l+l  l+l  l-l  l-l 


-  u    (p  n    -  p  n  ) 

i    k    l+l          i-l;   At        n/A+2^ 
"IS +   0(At    ) 


Hence,    it   can  be  seen  that 


M.        E.        -   M.    E.        ^___ 
l  l l      l       oPU 

At  +  d~F~  '  i 


.    Ax  +  Pn  En  U* 
ill 


P*     E*     U* 

l-l     l-l     l-l 


=  0(At2)   +  0(Ax2)   +   O(AtAx) 
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Therefore 


(pE)*+±  -    (pE)*        (pEtl)*  -    (pEU)*  (FU)*  -    (HI)  * 

-— — — —  + +  — — — - 

At  Ax  Ax 


dpE        SpEU       oTU        _/A    v        «/A.\        A/At    N 


(c)      OIL 

As  stated  previously,  the  mass  transport  terms  for  continuous 

~n  /  \ 

rezone  and  OIL  are  equal.  Also,  the  E.  (i  an  integer)  for  continuous 

rezone  and  OIL  are  equal  and  hence  the  calculations  of  the  order  of 

the  error  for  OIL  and  continuous  rezone  are  the  same. 
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I  ABSTRACT 


The  description  and  comparison  of  several  explicit  differencing 
schemes  for  the  solution  of  the  inviscid,  non-thermal-conducting,  single- 
material,  one-dimensional  equations  of  compressible  hydrodynamics  is  in- 
cluded in  this  report.   The  comparison  was  based  on  the  plots  obtained 
from  a  series  of  one  dimensional  test  problems  which  are  included  in 
APPENDIX  A. 
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